Establishment, FEM analysis and experimental validation of tooth movement prediction model of orthodontic archwire T-loop

Background The T-loop has been used clinically to close gap between teeth. And it is a typical orthodontic archwire bending method. However, the design of the T-loop parameters for different patients is based on the clinical experience of the dentists. The variation in dentists' clinical experience is the main reason for inadequate orthodontic treatment, even high incidence of postoperative complications. Methods Firstly, the tooth movement prediction model is established based on the analysis of the T-loop structure and the waxy model dynamic resistance. As well as the reverse reconstruction of the complete maxillary 3D model based on the patient CBCT images, the oral biomechanical FEM analysis is completed. A maxillary waxy dental model is manufactured to realize the water-bath measurement experiment in vitro mimicking the oral bio-environment. Thus, the calculated, simulation and experimental data are obtained, as well as obtaining a cloud of total deformation from the simulation analysis. Results The growth trend of the 11 sets of simulation data is the same as that of the experimental data. And all of them show that the tooth displacement is positively correlated with the cross-sectional size of the archwire, and the clearance distance. As well as the higher Young's modulus of the archwire material, the greater the tooth displacement. And the effect of archwire parameters on tooth displacement derived from simulation and experimental data is consistent with the prediction model. The experimental and calculated data are also compared and analyzed, and the two kinds of data are basically consistent in terms of growth trends and fluctuations, with deviation rates ranging from 2.17  to  10.00%. Conclusions This study shows that the accuracy and reliability of the tooth movement prediction model can be verified through the comparative analysis and deviation calculation of the obtained calculated, simulation and experimental data, which can assist dentists to safely and efficiently perform orthodontic treatment on patients. And the FEM analysis can achieve predictability of orthodontic treatment results.


Background
In recent years, oral malocclusion has shown a high prevalence and is one of the three major oral diseases [1]. It affects oral function, craniofacial development,

Open Access
*Correspondence: jiangjingang@hrbust.edu.cn appearance and interpersonal relationships, and even causes psychological disorders in minors [2,3,4]. Orthodontic techniques can effectively treat malocclusion. During orthodontic treatment, due to the phenomenon of orthodontic crowding, some teeth are selectively extracted and a corresponding extraction gap is formed. For this reason, the closure of the extraction gap is a very important aspect [5], which is directly related to the smooth implementation of the entire orthodontic plan [6,7,8].
T-loop is widely used to close gap between teeth due to the structural characteristics that can generate horizontal, vertical and torsional force [9,10,11,12,13]. During the process of closing gap, the teeth are subjected to both the orthodontic force generated by the T-loop and the biological resistance of the periodontal tissues [14,15,16]. In turn, the orthodontic force generated by the archwire are inextricably linked to the clearance distance, cross-sectional size and material properties of the T-loop [17,18].
In traditional orthodontic diagnosis, depending on the treatment stage and specific orthodontic purpose, dentists select archwire parameters empirically [19,20,21,22]. This leads to unpredictable tooth movement, and can easily cause irreparable damage to the patients and reduce treatment efficiency [13]. During the treatment, if the tooth movement is too small, the orthodontic effect is not obvious, and if the tooth movement is too large, it can increase the patient's dental pain, leading to periodontal tissue damage, tooth loosening, and even serious root and alveolar bone resorption problems [23,24,25]. Therefore, dentists will find the optimal position by repeatedly adjusting the parameters of the orthodontic archwire. However, repeated adjustments can reduce the mechanical properties of the orthodontic archwire. And it also increases the clinical treatment time and causes inconvenience to both the dentist and the patient. Therefore, by revealing the relationship between archwire parameters, orthodontic force and tooth movement, it enables the dentist to develop an individualized orthodontic treatment plan by taking into account the patient's malocclusion, the mechanical behavior of archwire and the tooth movement pattern. This further improves the efficiency and safety of fixed orthodontic treatment and reduces the patient's pain.
In this paper, based on the dynamic resistance model of the waxy model and the dynamic orthodontic force model of the T-loop, a tooth movement prediction model is established to realize the tooth movement prediction analysis. The T-loop parameters and tooth movement distance are quantified. And FEM analysis and experimental validation are both conducted. Experimental data are collected using the designed tooth movement measurement method. And the accuracy of the tooth movement prediction model is demonstrated by data comparison and deviation analysis. Therefore, the tooth movement prediction model established in this paper can assist dentists in selecting the appropriate T-loop parameters for the movement distance and angle of the tooth. Thus, it provides a theoretical basis for the development of individualized orthodontic treatment plans for patients. At the same time, it reduces the clinical treatment time of patients and dentists, improves orthodontic efficiency, as well as achieves predictability of orthodontic treatment effect.

Dynamic resistance analysis of waxy model
In order to restore the relative tooth movement under realistic conditions, a non-rigid waxy model is chosen to simulate the biological properties of the periodontal environment. Non-rigid waxy models are widely used in dental clinical experiments and training because the viscous fluid properties vary with temperature, such as the Typodon orthodontic training model [26,27]. As well as Liu et al. experimentally verified the feasibility of the waxy model for simulating tooth movement at a certain temperature [28]. Meanwhile, according to the theory of viscous fluid dynamics, during orthodontic treatment, the impedance force of tooth in waxy model is composed of resistance and inertia force [29]. Through experimental verification, the optimal water-bath temperature for tooth movement in the waxy model is determined to be 75 °C [30,31]. At the same time, an expression for the variation of the waxy dental density with time is established. According to the expression, the resistance equation and the inertia force equation, the dynamic resistance model of the waxy dental model in the process of simulated tooth movement is established [32,33], as shown in Fig. 1.

Modeling of dynamic orthodontic force in T-loop
During orthodontic treatment, the tooth movement is influenced by a combination of the periodontal tissue's biological resistance and the orthodontic force generated by the archwire. Besides, orthodontic force changes dynamically as treatment progresses [32]. Therefore, before establishing a tooth movement prediction model, a dynamic orthodontic force model needs to be established.
Analysis of T-loop working as follows: the dentist pulls the horizontal arm of T-loop to create a clearance, causes elastic deformation of the vertical arm and the arc portion, and then fixes the horizontal arm on the bracket. The tooth is pulled to a predetermined position using the restoring force of the T-loop itself, thus completing the closure of the gap between teeth [19,20]. At the same time, T-loop with symmetrical structure is used, that is, the Alpha/Beta structure on both sides of T-loop is equal, so that the orthodontic force released by the horizontal arm can be equal and opposite [34]. Therefore, only one side of T-loop needs to be analyzed. The T-loop deformation analysis is shown in Fig. 2.
The orthodontic force generated by the T-loop is composed of the vertical arm and the arc portion of the deformation superimposed. So the vertical arm and the arc portion are modeled mechanically separately. And then the total orthodontic force of the T-loop is calculated by superposition. Firstly, the vertical arm is subjected to force analysis and mechanical modeling. As shown in Fig. 2, the bending radius of the arc portion is R, the total height is h, the length of the vertical arm is y. By pulling the horizontal arm of the T-loop creates a clearance  The approximate differential equation for the vertical arm deflection curve is as follows: where E is the elastic modulus of material, M(x) is the bending moment of the vertical arm at x distance. I z is the rotational inertia of the archwire z-axis. The rotational inertia of a round archwire is I z = πD 4 /64, and D is the diameter of the round archwire. For a rectangular archwire I z = c 1 c 2 3 /12. The c 1 is length of the side perpendicular to the z-axis on the rectangular archwire cross-section. The rotation angle equation θ(x) and the deflection equation ν(x) for the vertical arm can be obtained by integrating Eq. (1) once and twice, respectively.
where C 0 and D 0 are integration constants that can be determined by the boundary conditions at special points. The expression for the bending moment of the vertical arm is: Substituting Eqs. (1) and (3) into Eq. (2) for integration, obtaining: Using the special points with known deflection and rotation angle, the integration constants in Eq. (4) can be determined. Since the longitudinal symmetry plane exists at the intersection of horizontal arm of the arc portion and vertical arm, and the symmetry plane is subject to external force, the axis of the deformed curved beam is still located within the longitudinal symmetry plane. The intersection point deformation belongs to the plane bending deformation of the curved beam. Therefore, the arc at the intersection can be equated to a curved beam with radius, curvature and section diameter. To take the radius micro-element dα cross section, as shown in Fig. 3.
After deformation, the differential equation of the deflection curve at the intersection of horizontal arm of the arc portion and vertical arm is as follows: where u is the displacement of the curved beam in the x-direction intersection cross section. The bending moment condition at the intersection is known to be M 0 = M| x=0 = −Py . I ω is the rotational inertia of the curved beam's cross section on the ω axis. The vertical arm is consistent with the bending form of the curved beam at the intersection, so I ω = I z . Also it is defined the arc length equation ds = Rdα, and obtained the differential equation of the deflection curve at the intersection point after substitution, as shown in Eq. (6): The non-associative differential equation with constant coefficients for the deflection curve after deformation of the curved beam at the intersection point is: The boundary condition at the intersection is u| α=π/2 = 0 , du dα α=π/2 = 0 . Therefore, it can be deter- EI ω , and bring them into Eq. (7) to obtain Eq. (8):

Fig. 3 Deformation analysis of nodes and horizontal arm of the arc portion
Also, the corner equation of the curved beam can be obtained as: (4) and (10), obtaining: According to the established deflection differential equation, the corresponding boundary conditions are replaced: the maximum deflection value m, and the maximum rotation angle at x = y. The external force that deforms the vertical arm is then obtained as: Further according to the principle of action and reaction force, the force P that deforms the vertical arm is equal to the counter force of the orthodontic force F 1 applied to the tooth during the orthodontic procedure. Thus, it can be obtained that: As mentioned earlier, the arc portion and the deformation of the vertical arm together constitute the orthodontic force generated by the T-loop during correction. After the analysis of the vertical arm, a mechanical modeling.
analysis of the arc portion is performed. Similar to the analysis of the vertical arm, the horizontal arm of the arc portion as the main deformation part, need to establish the corresponding deformation equation, as shown in Fig. 3.
During the deformation of the T-loop, the position of the connection between the horizontal arm of the arc portion and the vertical arm is constantly changing. It is not easy to measure the distance that the horizontal arm of the arc portion moves. Therefore, it is necessary to simplify the movement of the intersection symmetry center during the T-loop deformation. The symmetry center of the intersection before deformation is overlapped, and the difference between the theoretical length and the initial length of the vertical arm after deformation is calculated. This difference is the bending deformation s of the horizontal arm of the arc portion, which can be obtained from the following equation: where y is the length of the vertical arm in the natural state and also the bending deflection of the vertical arm. Similar to Eq. (1), the deflection curve differential equation for the horizontal arm of the arc portion is: The corner Eq. (15) and the deflection Eq. (16) for the horizontal arm of the arc portion can be determined by integrating Eq. (14) once and twice, respectively: where C 1 and D 1 are integration constants, and the bending moment expression for the horizontal arm of the arc portion is: where G is the component that deforms the horizontal arm of the arc portion, w is the total length of the horizontal arm of the arc section, and R is the bending radius of the arc portion. Substituting Eqs. (14) and (17) into Eqs. (15) and (16) for integration, the following equation is obtained: The lateral arc on one side can be equated to a curved beam of radian π/4, and the radian dβ is taken to be infinitesimal. This is the same process as solving for the vertical arm boundary condition, so the boundary condition can be given as: The boundary conditions of the arc portion are u| β=π = 0 , du dβ β=π = 0 , and the general solutions Then the deflection equation of the arc portion is: The corner equation of the arc portion is: Because the boundary conditions are v| l=0 = u| β=0 = 0 and θ | l=0 = ε| β=0 = 0 , it is able to obtain C 1 = 0 and D 1 = 0. Substituting C 1 = 0 and D 1 = 0 into Eqs. (18) and (19) obtains Eqs. (23) and (24): The corresponding boundary conditions are substituted as follows: when l = w-R, the horizontal arm of the arc portion reaches the maximum deflection at its end, which is v = s, approximated by the same Eq. (13). Further can get the bending force of the horizontal arm of the arc portion, as follows: From further mechanical analysis, it is obtained that the component force along the horizontal arm should be equal to the orthodontic force F 2 released by the arc portion. The following equation is obtained: From the mechanical analysis, the restoring force from the deformation of the vertical arm and the arc portion together constitute the orthodontic force F 0 released by the horizontal arm. That is, the following equation:

Establishment of tooth movement prediction model
After the previous calculations, the dynamic resistance model for the tooth movement simulated in the waxy model can be expressed by Eq. (28): where f is the dynamic resistance to tooth movement in the waxy model. Bringing in the orthodontic force F 0 released by the T-loop, it is obtained that: In physics, the relationship between the displacement r and the velocity vector v, Eq. (30), can be determined by integrating: In Newton's Second Law, the relationship between the acceleration a 0 , the mass m 0 and the force F can be expressed by Eq. (32): The acceleration a 0 of the tooth moving in the waxy model, subject to orthodontic force F 0 , inertial force f L and resistance f D , can be expressed by Eq. (32): By substituting Eq. (33) into (31), the tooth movement prediction model can be obtained. The following equation can be obtained: where S m is the amount of tooth movement under the action of the T-loop in the waxy model, m 0 is the mass of the tooth, and t 1 is the time of tooth movement.

Tooth movement FEM analysis
In this study, a 3D model of the complete maxillary bone, dental row and periodontal ligament is reconstructed in reverse by using CBCT images (EWOO-VATECH Co., Ltd., Korea, scan thickness 0.5 mm, 301 scanned images) based on the patient's skull. And according to the clinical orthodontic treatment process, a 3D model with brackets (28) and T-loop is installed on the reconstructed dental. At last, the orthodontic biomechanical FEM analysis is performed for the complete assembly model. The flowchart of the reverse 3D reconstruction model is shown in Fig. 4. And the simulation cloud is shown in Fig. 5. In the simulation process, the archwire is adopted of stainless steel archwire and Australian archwire, which are commonly used in clinical orthodontics. In addition, the material properties of the brackets and maxillary components are shown in Table 1 [34,35,36].
The simulation parameters shown in Fig. 5 are: (1) The archwire is a 0.014-inch round Australian archwire. respectively. In addition, the displacement load is defined as follows: set 10 analysis steps so that the clearance distance grows linearly with 0.1 of the step length, and take the clearance distance as 0.3-1.2 mm for FEM analysis. And the maximum displacement of NO. 13 tooth is 1.4575 mm at the cusp and the minimum displacement is 0.2512 mm at the root, which meets the clinical orthodontic requirements [37,38].

Tooth movement experimental measurement method
In this study, the waxy model is used to simulate the periodontal tissue structure. Firstly, alginate is used as the impression material, and the separated teeth are placed into the alginate negative mold, then the molten waxy base is poured, left to cool and removed. Thus, the preparation of waxy dental model for experimental use is completed. Due to the viscoelastic properties at the corresponding temperature, the waxy model under the action of external force will produce a corresponding elastic deformation [28]. Finally, the teeth on the waxy model will complete the simulation of the orthodontic process with the change of time. The experimental system composition and the waxy dental model are shown in Fig. 6. Also detailed waxy model preparation process and modeling methods have been described in the study of Huang et al. [32]. The position of each maxillary teeth is marked using the FDI (Fédération Dentaire Internationale) marking method, and maxillary teeth No. 13 and No. 15 are selected as the teeth for the experimental study. Before clinical orthodontic treatment, a total of four teeth in the maxilla and mandible are usually extracted to allow enough space for remaining teeth to move. In the experiment, tooth No.14 is extracted and the two horizontal arms of the T-loop are ligated into the brackets of teeth No. 13 and No. 15. The highest points of the marginal crests of the two teeth are marked, as shown in Fig. 6. The distance between the two marker points is measured before the water-bath experiment, and then the waxy model is immersed in a constant temperature water-bath at 75 °C for 2 min. After taking out, the distance between the two marker points is measured No.13

Fig. 6
Installed waxy dental model with T-loop and tooth movement experiment system again with a vernier caliper [32]. The difference of the distance between the two marker points before and after the water-bath is sought as the experimental value of the tooth movement distance under the action of T-loop.

Results
In clinical treatment, the magnitude of orthodontic force is changed by adjusting the characteristic parameters of the T-loop. And the magnitude of orthodontic force has an important influence on the tooth movement, the periodontium disease, tooth root resorption, and the alveolar bone reconstruction [39]. Therefore, in this study, the accuracy and reliability of the tooth movement prediction model is verified by comparing and analyzing the calculated, simulation and experimental data. In the FEM analysis, 11 kinds of T-loop 3D models with different parameters are established based on clinical orthodontic treatment. The loading process of the FEM analysis is the same as described before, with the aim of investigating the effects of different clearance distances, cross-sectional and materials on tooth movement. To be easily represented, the orthodontic archwires are codenamed as follows: I-S1616, II-S1622, III-S1625, IV-S0014, V-S0016, VI-S0018, VII-S0020, VIII-A0014, IX-A0016, X-A0018, XI-A0020. Where S is stainless steel archwire, A is Australian archwire, 1616 indicates 0.016 × 0.016 inch rectangular archwire, and 0014 indicates 0.014 inch round archwire. Similarly, the 11 kinds of T-loop bending parameters used in the experimental are the same as those of the 3D model, and 10 measurement points are set at 0.1 mm intervals to collect data during the loading measurement stage of the experiment. The calculated, simulation and experimental data are shown in Table 2.
As well, the total strain clouds for the biomechanical FEM analysis of tooth movement under the action of 11 kinds of T-loop are shown in Fig. 7.
In this paper, the simulation data and experimental data are compared and analyzed, so that the accuracy and reliability of simulation analysis and experimental measurement can be verified with each other. Figure 8 shows the effect of rectangular archwire with different cross-sectional widths (Codenamed as I, II, and III) on the variation of tooth displacement with T-loop clearance distance. Both simulation and experimental results shows that the tooth displacement is positively correlated with the cross-sectional width for the same clearance distance. Similarly, the tooth displacement is positively correlated with the clearance distance for the same cross-sectional width. In addition, the effect of cross-sectional width on tooth displacement increases as the clearance distance increases.
Comparing different round cross-sectional diameter (Codenamed as IV, V, VI, and VII), the effect of tooth displacement with the variation of T-loop clearance distance is shown in Fig. 9. Both simulation and experimental results show that the tooth displacement is positively correlated with the cross-sectional diameter of the T-loop for the same clearance distance, which is similar to the conclusion for rectangular archwire. From the simulation and experimental results, although the cross-sectional diameter of the archwire affects the tooth displacement, its effect is weaker than that of the clearance distance of the T-loop.
In clinical orthodontic treatment, stainless steel archwire and Australian archwire are mainly used. Therefore, two sets of results based on material are divided. The data of four identical archwire types in the two sets are compared to investigate the effect of different archwire materials on tooth displacement [5,32]. The stainless steel archwires are codenamed as IV, V, VI, and VII, and the Australian archwires are codenamed as VIII, IX, X, and XI. As shown in Fig. 10, both the simulation and experimental data show that the effect of tooth displacement with the same cross-sectional dimensions is more pronounced with increasing clearance distance using stainless steel archwire than Australian archwire. In addition, this conclusion became more evident as the diameter of the archwire increased.
While mutually verifying the accuracy and reliability of the simulation analysis and experimental measurements, this paper analyses the calculated and experimental data, and performs deviation calculations to verify the accuracy and validity of the tooth movement prediction model. The calculated data is shown in Table 2. And in Figs. 11, 12, 13, the experimental and calculated data, and deviation rates are shown by solid, dashed and dotted lines, respectively.
A comparison of experimental and calculated data of tooth displacement versus clearance distance for rectangular stainless steel archwire of different cross-sectional widths (Codenamed as I, II and III) is shown in Fig. 11. The obtained calculated and experimental data have the same trend and the deviation varies from 2.17 to 8.86%. When other conditions are the same, the trend of the curves shows that the calculated data of tooth displacement are positively correlated with the cross-sectional width, and the effect of the cross-sectional width on tooth displacement increases with the increase of the clearance distance.
For the round archwire, the comparison of experimental and calculated data of tooth displacement for different cross-sectional sizes of stainless steel archwire (Codenmed as IV, V, VI, VII) and Australian archwire (codenamed as VIII, IX, X, XI) at different clearance distances, and deviation rates are shown in Figs. 12 and 13. The trends of the calculated and experimental data are consistent in different archwire materials. The deviation rate of calculated and experimental data for stainless steel archwire ranges from 2.36 to 10.00%, while the deviation rate of Australian archwire ranges from 2.47 to 9.38%. For the two kinds of orthodontic archwire with different materials, the larger the cross-sectional size of the archwire, the greater the tooth displacement effect produces when the clearance distance is the same. When the crosssectional size are all the same, the tooth displacement effect produced by the round orthodontic archwire is positively correlated with the clearance distance.

Discussion
In order to reveal the tooth movement characteristics under different force during clinical orthodontic procedure, researchers have studied it through experimental measurements, FEM analysis and theoretical derivation. For example, researchers have established FEMs of tooth movement under different orthodontic force [40,41,42,43]. During the simulation, the magnitude and direction of the orthodontic force is artificially set without reference, which is not consistent with the clinical reality. In addition, periodontal tissues cannot be effectively differentiated by materials, so corresponding experimental validation should be performed. However,  experimental materials are mostly rigid models, such as plaster models or metal models [44], which lack biological properties. In addition, the orthodontic process is a dynamic one, and the orthodontic force will become weaker with the reduction of archwire deformation. Therefore, the real orthodontic effect cannot be effectively reproduced and the prediction function is not yet perfect. In contrast, the T-loop is a common orthodontic arch for clinical treatment of such oral malocclusion as closing the gap between teeth. Although the relationship between archwire shape and orthodontic force for T-loop has been studied to some extent [45,46,47,48], but the archwire bending parameter and orthodontic force relationship, orthodontic force and tooth movement relationship have not been established. Therefore, the current studies still cannot assist dentists to effectively develop individualized T-loop parameter designs for different patients, as well as to accurately predict orthodontic treatment effect. Thus, this paper firstly establishes a tooth movement prediction model based on the T-loop analysis and the waxy model dynamic resistance. Subsequently, to verify the accuracy of the prediction model, a reverse reconstruction model of the maxillary and dentition based on the patient's CBCT images is used to achieve a biomechanical FEM analysis of the tooth movement. Finally, a waxy dental model is established, and a water-bath measurement experiment mimicking the biological environment of the oral cavity is realized. Simulation analysis and experimental measurements are compared to verify each other's accuracy and reliability. And as shown in Figs. 8 and 9, the conclusion about the influence of different archwire parameters on tooth movement derived from the simulation analysis and experimental measurements is consistent with the established tooth movement prediction model. This further illustrates the accuracy of the simulation analysis and experimental measurements.
Afterwards, the prediction model and experimental measurements are compared and analyzed, and the deviation is calculated, and the total deviation range between the calculated and experimental data is obtained: 2.17 ~ 10.00%, which is within the acceptable range and is approved by the dentists of Peking University School of Dentistry. In addition, it can be determined from the trends shown in Figs. 11, 12 and 13 that the movement of the target tooth is positively correlated with both square and round archwire in terms of clearance distance and cross-sectional parameters. And the tooth displacement in the stainless steel archwire move faster than those in the Australian archwire as the clearance distance increases. The tooth displacement of the Australian archwire is milder for the same T-loop parameters. These trends and correlations are corrected with the existing academic studies in this field.

Conclusions
In this paper, based on the analysis of the T-loop structure and the waxy model dynamic resistance, a tooth movement prediction model is established. The relationship between T-loop parameters, orthodontic force and tooth movement are revealed and quantified. The validity and reliability of the simulation analysis and experimental measurements, as well as the accuracy of the tooth movement prediction model, are verified by comparing and analyzing the simulation analysis and experimental measurements with each other, as well as calculating and analyzing the deviation of the experimental and calculated data.  Therefore, it can be affirmed that the established tooth movement prediction model has a high degree of confidence and can assist dentists in selecting the appropriate archwire parameters when bending T-loop to obtain the desired tooth movement distance. At the same time, the tooth movement prediction model provides dentists with an effective theoretical guide to avoid inefficient bending of archwire and improve the treatment effect of malocclusion. Moreover, the process and results of the FEM analysis enable predictive orthodontic treatment effects, thus laying the foundation for the realization of digital orthodontic treatment.
In future work, other detailed factors that may affect tooth movement during orthodontic treatment will be further explored. Such as the relative static friction and relative sliding friction between the orthodontic archwire and brackets, and how to choose the right archwire parameters to achieve the best orthodontic effects and speed up the treatment process.  Comparison of experimental and calculated data with different cross-sectional widths, and deviation rate Fig. 12 Comparison of experimental and calculated data of stainless steel archwire at different clearance distances, and deviation rate